/* Evaluate recovered beliefs from GSU beliefs data on Raven) */

* log file
capture log close _all
log using "Evaluate Recovered Beliefs on Raven.log", replace name(recover)

* configure Stata
capture: version 18
set processors 4
set more off
set scheme s1color
set seed 987654321
timer clear 1
timer on 1

* tasks
global doIND "y"
global doPOOL "y"

* graphics font
graph set window fontface "Garamond"
graph set window fontface "Candara"

* create directories needed
capture: mkdir figures

* tell us what version ran
about

* define the symbol used for the CRRA utility function (not always r)
local crra "s"
local crra "r"

* number of subjects
local Nsub = 407

* name of file stem
local belief "raven"

* evaluate over subjects and questions
if "$doIND" == "y" {

qui {

noi: di "Evaluating beliefs of `Nsub' subjects..."

* loop over all subjects
forvalues s=1/`Nsub' {

* loop over all questions
forvalues q=1/36 {

local bin = 8

* get the recovered beliefs
use recovered/`belief'_s`s'_q`q', clear
summ

* label for mcmc sample
rename mcmc mcmc_size
generate long mcmc = _n
label variable mcmc "MCMC sample number"

* check beliefs sum to 1
generate check_belief = 0
forvalues x=1/`bin' {
	replace check_belief = check_belief + b_`x'
}
summ check_belief
rename check_belief normalize_b
label variable normalize_b "Normalization for recovered beliefs to sum to 1"

* get posterior predictive mean and median for each bin
generate sum_p_mean = 0
generate sum_p_median = 0
forvalues x=1/`bin' {
	summ b_`x', detail
	local m = r(mean)
	generate p_mean_`x' = `m'
	label variable p_mean_`x' "Posterior predictive mean for bin `x'"
	local m = r(p50)
	generate p_median_`x' = `m'
	label variable p_median_`x' "Posterior predictive median for bin `x'"
	replace sum_p_mean = sum_p_mean + p_mean_`x'
	replace sum_p_median = sum_p_median + p_median_`x'

	label variable r_`x' "Report for bin `x'"
	label variable b_`x' "Recovered belief for bin `x'"
}

* normalize posterior mean and median
forvalues x=1/`bin' {
	replace p_mean_`x' = p_mean_`x'/sum_p_mean
	replace p_median_`x' = p_median_`x'/sum_p_median
}
summ p_mean_*
summ p_median_*

rename sum_p_mean normalize_mean
label variable normalize_mean "Normalization for recovered posterior predictive mean to sum to 1"
rename sum_p_median normalize_median
label variable normalize_median "Normalization for recovered posterior predictive median to sum to 1"

* illustrate
des mcmc r_* b_* normalize_b normalize_mean p_mean_* normalize_median p_median_*
list mcmc r_* r eta phi in 1/5, noobs
list mcmc b_* normalize_b in 1/5, noobs
list mcmc normalize_mean p_mean_* in 1/5, noobs
list mcmc normalize_median p_median_* in 1/5, noobs

* summarize risk preference parameters and display
foreach x in r eta phi {
	summ `x'
	local `x' = r(mean)
}

* maximum likelihood parameters and beliefs
summ ll
local maxLL = r(max)
di "Maximum LL is `maxLL'"
summ ll, detail

* get the second highest LL value, and use it to get the maxLL (avoids rounding error)
sort ll
summ ll
local nobs = r(N)
list ll if _n == `nobs'-1
list ll if _n == `nobs'
summ ll if _n == `nobs'-1
local maxLL = r(mean)
di "Maximum LL is strictly greater then `maxLL', and here it is..."
summ ll if ll > `maxLL'

* now get the ML estimates
list r eta phi ll if ll == `maxLL', noobs
foreach x in r eta phi {
	summ `x' if ll > `maxLL'
	local `x'_ml = r(mean)
	di "ML estimate of `x' is ``x'_ml'"
}
forvalues x=1/`bin' {
	 summ b_`x' if ll > `maxLL'
	 generate ml_`x' = r(mean)
}
summ ml_*

* save, reduce to one observation and save the summary
summ
save recovered/`belief'_s`s'_q`q'_save, replace
keep if _n == 1

* get deviation of beliefs from reports
generate dev_p_mean = 0
forvalues x=1/`bin' {
	replace dev_p_mean = dev_p_mean + abs(r_`x' - p_mean_`x')
}
replace dev_p_mean = dev_p_mean*100
label variable dev_p_mean "Percent deviation of posterior mean beliefs from reports"
summ dev_p_mean
local dev = r(mean)
local dev_ = string(`dev', "%2.0f") + "%"
list r_* p_mean_* dev_p_mean, noobs
summ
save recovered/`belief'_s`s'_q`q'_summary, replace

* close loop over questions
}

* close loop over subjects
}

* end of qui
}

* end of $doIND
}


* collate files
if "$doPOOL" == "y" {

* qui block
qui {

noi: di "Collating individual files..."

* keep track of files
local f = 0

forvalues s=1/`Nsub' {
forvalues q=1/36 {

local bin = 8

* get the recovered and processed beliefs -- summary file
use recovered/`belief'_s`s'_q`q'_summary, clear
local f = `f'+1
save tmp`f', replace

* get the recovered and processed beliefs -- detailed save file
use recovered/`belief'_s`s'_q`q'_save, clear
keep id question r eta phi mu mcmc r_* b_* uid
save tmp`f'_, replace

}
}

noi: di "Appending individual files..."

* append summary file
use tmp1, clear
forvalues x=2/`f' {
	append using tmp`x'
	erase tmp`x'.dta
}
erase tmp1.dta
save `belief'_summary, replace

* append save file
use tmp1_, clear
forvalues x=2/`f' {
	append using tmp`x'_
	erase tmp`x'_.dta
}
erase tmp1_.dta
save `belief'_save, replace

* end of qui
}

* review
tab id
tab question

* end of $doPOOL
}



* time taken overall
timer off 1
timer list
local secs = r(t1)
local mins = `secs'/60
local hrs = `mins'/60
local secs_ = string(`secs', "%10.0f")
local mins_ = string(`mins', "%4.1f")
local hrs_ = string(`hrs', "%4.2f")
di "Complete calculations took `secs_' seconds, `mins_' minutes, or `hrs_' hours."

log close recover
